# R version 4.5.0
# Run '1. dataset' first.

library(haven)
library(dplyr)
library(rcompanion)
library(tidyr)
library(ggplot2)
library(cowplot)
library(reshape2)
library(scales)
library(extrafont)
library(ggtext)
library(MASS)
library(sfsmisc)
library(mediation)
library(patchwork)

setwd() # set working directory

# Figure2-3
# Download raw data from: https://kgss.skku.edu/kgss/data.do
kgss <- read_sav("★2003-2023_KGSS_kor_public_v3.sav")
kgss <- kgss %>% filter(!PARTYLR == -8)

kgss$IDEO <- ifelse(kgss$PARTYLR %in% c(1,2), "L", 
                    ifelse(kgss$PARTYLR %in% c(4,5), "R", "I"))

for (year in 2010:2023) {
  assign(paste0("kgss", substr(year, 3, 4)), 
         kgss %>% filter(YEAR == year))
}

kgss_list <- list(kgss10, kgss11, kgss12, kgss13, kgss14, kgss16, kgss18, kgss21, kgss23)
years <- c("10", "11", "12", "13", "14", "16", "18", "21", "23")
filter_prtyid <- function(df, year) {
  prtyid_var <- paste0("PRTYID", year)
  df %>% filter(!.data[[prtyid_var]] %in% c(-8, 77, 88, 99))
}
kgss_list <- mapply(filter_prtyid, kgss_list, years, SIMPLIFY = FALSE)

list2env(setNames(kgss_list, paste0("kgss", years)), envir = .GlobalEnv)

cramers_v_results <- list()
cramers_v_results2 <- list()

for (i in seq_along(kgss_list)) {
  year <- years[i]
  df <- kgss_list[[i]]
  
  prtyid_var <- paste0("PRTYID", year)
  
  contingency_table <- table(df[[prtyid_var]], df$PARTYLR)
  contingency_table2 <- table(df[[prtyid_var]], df$IDEO)
  cramers_v <- cramerV(contingency_table)
  cramers_v2 <- cramerV(contingency_table2)
  
  cramers_v_results[[year]] <- cramers_v
  cramers_v_results2[[year]] <- cramers_v2
}

cramers_v_df1 <- data.frame(
  Year = as.numeric(names(cramers_v_results)) + 2000,
  CramerV = sapply(cramers_v_results, as.numeric),
  Dataset = "Dataset1"
)

cramers_v_df2 <- data.frame(
  Year = as.numeric(names(cramers_v_results2)) + 2000,
  CramerV = sapply(cramers_v_results2, as.numeric),
  Dataset = "Dataset2"
)

cramers_v_df <- rbind(cramers_v_df1, cramers_v_df2)

f3b <- ggplot(cramers_v_df, aes(x = Year, y = CramerV, color = Dataset, group = Dataset)) +
  geom_line(linewidth = 1.05) +
  geom_point() +
  labs(title = "<b>Partisan Sorting Measured by Cramer's V<b>",
       x = " ",
       y = " ",
       color = "Type",
       caption = "") + 
  theme_minimal() +
  theme(
    plot.title = element_markdown(size = 14, lineheight = 1.1, family = "roboto"),
    plot.margin = margin(t = 20, r = 10, b = 10, l = 10), 
    legend.position = "bottom" 
  ) +
  scale_color_manual(values = c("Dataset1" = "#c0c0c0", "Dataset2" = "#8a6e96"),
                     labels = c("Ideology in 5 scale", "Ideology in 3 scale"))  # Custom legend labels

addmargins(table(kgss10$PRTYID10, kgss10$IDEO)); p10d <- 185/1063; p10p <- 229/1063

addmargins(table(kgss11$PRTYID11, kgss11$IDEO)); p11d <- 147/977; p11p <- 225/977

p12d <- addmargins(table(kgss12$PRTYID12, kgss12$IDEO))[2,2]/addmargins(table(kgss12$PRTYID12, kgss12$IDEO))[9,4]; p12p <- addmargins(table(kgss12$PRTYID12, kgss12$IDEO))[1,3]/addmargins(table(kgss12$PRTYID12, kgss12$IDEO))[9,4]

addmargins(table(kgss13$PRTYID13, kgss13$IDEO)); p13d <- 165/890; p13p <- 274/890

addmargins(table(kgss14$PRTYID14, kgss14$IDEO)); p14d <- 188/840; p14p <- 252/840

addmargins(table(kgss16$PRTYID16, kgss16$IDEO)); p16d <- 130/767; p16p <- 166/767

addmargins(table(kgss18$PRTYID18, kgss18$IDEO)); p18d <- 299/771; p18p <- 91/771

addmargins(table(kgss21$PRTYID21, kgss21$IDEO)); p21d <- 214/772; p21p <- 203/772

addmargins(table(kgss23$PRTYID23, kgss23$IDEO)); p23d <- 176/793; p23p <- 311/793

data <- data.frame(
  Year = c(2010, 2011, 2012, 2013, 2014, 2016, 2018, 2021, 2023),
  p_d = c(p10d, p11d, p12d, p13d, p14d, p16d, p18d, p21d, p23d),
  p_p = c( p10p, p11p, p12p, p13p, p14p, p16p, p18p, p21p, p23p)
)

data$p_rest <- 1 - data$p_d - data$p_p

data_long <- data %>%
  pivot_longer(cols = c(p_d, p_p, p_rest), names_to = "Category", values_to = "Proportion")

data_long$Category <- factor(data_long$Category, levels = c("p_rest", "p_d", "p_p"))

custom_colors <- c("p_p" = "#7851A9", "p_d" = "#DDA0DD", "p_rest" = "grey80")

f3a <- ggplot(data_long, aes(x = factor(Year), y = Proportion, fill = Category)) +
  geom_bar(stat = "identity") +
  labs(
    title = "<b>Proportion of Sorted and Unsorted Partisans<b>",
    x = " ",
    y = " ",
    caption = "Source: Korean General Social Survey (2010-2023)",
    fill = "Type"
  ) +
  scale_y_continuous(labels = scales::percent) + 
  scale_fill_manual(
    values = custom_colors, 
    labels = c("p_p" = "Conservative-sorted", "p_d" = "Democrat-sorted", "p_rest" = "Unsorted partisans")
  ) +
  theme_minimal() +
  theme(
    plot.title = element_markdown(size = 14, lineheight = 1.1, family = "Roboto"),
    plot.margin = margin(t = 20, r = 10, b = 10, l = 10), 
    plot.caption = element_text(size = 8, color = "grey50", hjust = 0),
    panel.grid = element_blank(),
    legend.position = "bottom" 
  )

f3a + f3b
ggsave("figure3.jpeg", width = 10, height = 6, dpi = 300)

proportion_LR <- list()

for (i in seq_along(kgss_list)) {
  year <- years[i]
  df <- kgss_list[[i]]
  
  prtyid_var <- paste0("PRTYID", year)
  
  contingency_table <- table(df[[prtyid_var]], df$IDEO)
  
  proportion_table <- prop.table(contingency_table, margin = 2)
  
  proportion_LR[[year]] <- proportion_table
}

melted_list <- list()

for (year in names(proportion_LR)) {
  melted_df <- melt(proportion_LR[[year]], variable.name = "Ideology", value.name = "Proportion")
  melted_df$year <- as.numeric(year)
  melted_df$rowname <- rownames(proportion_LR[[year]])
    melted_list[[year]] <- melted_df
}

proportion_df <- bind_rows(melted_list)
proportion_df$year <- proportion_df$year + 2000

proportion_df$party_affiliation <- case_when(
  proportion_df$rowname == "1" & proportion_df$year %in% c(2018, 2021, 2023) ~ "DEM",
  proportion_df$rowname == "2" & !(proportion_df$year %in% c( 2018, 2021, 2023)) ~ "DEM",
  proportion_df$rowname == "1" & proportion_df$year %in% c(2010, 2011, 2012, 2013, 2014, 2016) ~ "PPP",
  proportion_df$rowname == "2" & proportion_df$year %in% c(2018, 2021, 2023) ~ "PPP",
  TRUE ~ "Others"
)

proportion_L_df <- proportion_df %>% filter(Var2 == "L")
proportion_L_df$party_affiliation <- factor(proportion_L_df$party_affiliation, levels = c("Others", "PPP", "DEM"))
proportion_R_df <- proportion_df %>% filter(Var2 == "R")
proportion_R_df$party_affiliation <- factor(proportion_R_df$party_affiliation, levels = c("Others", "DEM", "PPP"))

plot_left <- ggplot(proportion_L_df, aes(x = factor(year), y = Proportion, fill = party_affiliation)) +
  geom_bar(stat = "identity", position = "fill") +
  labs(
    title = "Party Identification Among Liberals",
    x = " ",
    y = " ",
    fill = "Party Identification",
    caption = "Source: Korean General Social Survey (2010-2023)"
  ) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0, 1)) +
  scale_fill_manual(
    values = c("royalblue", "#CD5C5C", "lightgrey"),
    labels = c("Democratic parties", "Conservative parties", "Other parties"),
    breaks = c("DEM", "PPP", "Others")
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.margin = margin(t = 20, r = 10, b = 20, l = 10),
    plot.caption = element_text(hjust = 0, size = 8),
    legend.position = "none"  # Remove legend for the first plot
  )

plot_right <- ggplot(proportion_R_df, aes(x = factor(year), y = Proportion, fill = party_affiliation)) +
  geom_bar(stat = "identity", position = "fill") +
  labs(
    title = "Party Identification Among Conservatives",
    x = " ",
    y = " ",
    fill = "Party Identification",
    caption = "Source: Korean General Social Survey (2010-2023)"
  ) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0, 1)) +
  scale_fill_manual(
    values = c("royalblue", "#CD5C5C", "lightgrey"),
    labels = c("Democratic parties", "Conservative parties", "Other parties"),
    breaks = c("DEM", "PPP", "Others")
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.margin = margin(t = 20, r = 10, b = 20, l = 10),
    plot.caption = element_text(hjust = 0, size = 8)
  )

plot_left <- plot_left +
  labs(title = "Party Identification Among Liberals") +
  theme(
    plot.title = element_text(size = 25, face = "bold", hjust = 0.5),
    axis.title.y = element_text(size = 16),
    axis.text.y = element_text(size = 14),
    axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
    plot.caption = element_text(size = 12, hjust = 0),
    plot.margin = margin(t = 10, r = 20, b = 10, l = 20) # Reduce side margins
  )

plot_right <- plot_right +
  labs(title = "Party Identification Among Conservatives") +
  theme(
    plot.title = element_text(size = 25, face = "bold", hjust = 0.5),
    axis.title.y = element_text(size = 16),
    axis.text.y = element_text(size = 14),
    axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
    plot.caption = element_blank(),
    plot.margin = margin(t = 10, r = 20, b = 10, l = 20) # Reduce side margins
  )

legend_plot <- ggplot(proportion_R_df, aes(x = factor(year), y = Proportion, fill = party_affiliation)) +
  geom_bar(stat = "identity", position = "fill") +
  scale_fill_manual(
    values = c("royalblue", "#CD5C5C", "lightgrey"),
    labels = c("Democratic parties", "Conservative parties", "Other parties"),
    breaks = c("DEM", "PPP", "Others"),
    name = "Party Identification"  # Correct legend title
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.title = element_text(size = 17),  # Adjust legend title size
    legend.text = element_text(size = 14)   # Adjust legend labels size
  )

all_components <- cowplot::get_plot_component(legend_plot, "guide-box", return_all = TRUE)
legend <- all_components[[3]]  # Select correct legend

combined_plot <- plot_grid(
  plot_left + theme(legend.position = "none"),
  plot_right + theme(legend.position = "none"),
  ncol = 2, align = "hv", rel_widths = c(1, 1)
)

final_plot <- plot_grid(
  combined_plot, legend, ncol = 1, rel_heights = c(1, 0.05)  # Adjust height for legend
) +  theme(plot.margin = margin(t = 20, r = 10, b = 10, l = 10)) 

ggsave("figure2.jpeg", plot = final_plot, width = 17, height = 10, dpi = 300)

# Figure 4
# Originial data downloaded from: https://v-dem.net/data/the-v-dem-dataset/
ldi <- read.csv("Vdem_LDI_short.csv")

ldi_long <- ldi %>%
  pivot_longer(cols = -Year, names_to = "Country", values_to = "Value")

ldi_long <- ldi_long %>%
  mutate(Country = ifelse(Country == "Korea", "South Korea", Country))

ldi_plot <- ggplot(ldi_long, aes(x = Year, y = Value, linetype = Country, 
                               size = ifelse(Country == "South Korea", 1.8, 1.2))) +
  geom_line() +
  scale_size_identity() +
  scale_y_continuous(limits = c(0.3, 1)) +
  scale_x_continuous(breaks = c(2000, 2005, 2010, 2015, 2020, 2024)) +
  labs(title = "V-Dem Liberal Democracy Index Across Selected Countries") +
  xlab(NULL) + ylab(NULL) +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 13, hjust = 0.5, margin = margin(b = 10)),
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 10)
  )

print(ldi_plot)
ggsave("figure4.jpeg", plot = ldi_plot, width = 6.3, height = 4.4, units = "in", dpi = 300)

# Figure 5
apr20 <- read_sav("apr2020_df.sav") 
apr20.DEM <- subset(apr20, PID == "DEM")
apr20.MTD <- subset(apr20, PID == "MTD") # MTD = PPP, the biggest conservative party; it changed its name

set.seed(2024)
apr20.DEM.med.y <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.out.y <- lm(Y ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.med.out.y <- mediate(apr20.DEM.med.y, apr20.DEM.out.y, 
                               treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)
apr20.MTD.med.y <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.out.y <- lm(Y ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.med.out.y <- mediate(apr20.MTD.med.y, apr20.MTD.out.y, 
                               treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

png("figure5.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr20.DEM.med.out.y, main = "Democratic Party (April 2020)", cex.main = 0.95)
plot(apr20.MTD.med.out.y, main = "Conservative Party (April 2020)", cex.main = 0.95)
dev.off()

# Figure 6
apr22 <- read_sav("apr2022_df.sav") 
apr22.DEM <- subset(apr22, PID =="DEM")
apr22.PPP <- subset(apr22, PID =="PPP")

apr22.DEM.med.y <- lm(AP ~ SORTING + GENDER + AGE + EDUC + INCOME, data = apr22.DEM)
apr22.DEM.out.y <- lm(Y ~ AP + SORTING + GENDER + AGE + EDUC + INCOME, data = apr22.DEM)
apr22.DEM.med.out.y <- mediate(apr22.DEM.med.y, apr22.DEM.out.y, 
                              treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)
apr22.PPP.med.y <- lm(AP ~ SORTING + GENDER + AGE + EDUC + INCOME, data = apr22.PPP)
apr22.PPP.out.y <- lm(Y ~ AP + SORTING + GENDER + AGE + EDUC + INCOME, data = apr22.PPP)
apr22.PPP.med.out.y <- mediate(apr22.PPP.med.y, apr22.PPP.out.y, 
                              treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)
png("figure6.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr22.DEM.med.out.y,  main="Democratic Party (April 2022)", cex.main=0.95)
plot(apr22.PPP.med.out.y,  main="Conservative Party (April 2022)", cex.main=0.95)
dev.off()

# Figure 7
jun22 <- read_sav("jun2022_df.sav") 
jun22.DEM <- subset(jun22, PID =="DEM")
jun22.PPP <- subset(jun22, PID =="PPP")

jun22.DEM.med.y <- lm(AP ~ SORTING + AGE + GENDER + INCOME + EDUC, data=jun22.DEM)
jun22.DEM.out.y <- lm(Y ~ SORTING + AP + AGE + GENDER + INCOME + EDUC, data=jun22.DEM)
jun22.DEM.med.out.y <- mediate(jun22.DEM.med.y, jun22.DEM.out.y, sims=500, treat = "SORTING", mediator = "AP")
jun22.PPP.med.y <- lm(AP ~ SORTING + AGE + GENDER + INCOME + EDUC , data=jun22.PPP)
jun22.PPP.out.y <- lm(Y ~ SORTING + AP + AGE + GENDER + INCOME + EDUC , data=jun22.PPP)
jun22.PPP.med.out.y <- mediate(jun22.PPP.med.y, jun22.PPP.out.y, sims=500, treat = "SORTING", mediator = "AP")

png("figure7.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(jun22.DEM.med.out.y, main="Democratic Party (June 2022)",cex.main=0.95)
plot(jun22.PPP.med.out.y, main="Conservative Party (June 2022)",cex.main=0.95)
dev.off()

# Figure 8
apr24 <- read_sav("apr2024_df.sav") 
apr24.DEM <- subset(apr24, PID =="DEM")
apr24.PPP <- subset(apr24, PID =="PPP")

apr24.DEM.med.y <- lm(AP2~SORTING + AGE + GENDER+ INCOME + EDUC, data=apr24.DEM)
apr24.DEM.out.y <- lm(Y~SORTING + AGE + GENDER+ INCOME + EDUC + AP2, data=apr24.DEM)
apr24.DEM.med.out.y <- mediate(apr24.DEM.med.y, apr24.DEM.out.y, 
                          treat = "SORTING", mediator = "AP2", robustSE = T, sims = 500)
apr24.PPP.med.y <- lm(AP2~SORTING + AGE + GENDER + INCOME + EDUC, data=apr24.PPP)
apr24.PPP.out.y <- lm(Y~SORTING + AGE + GENDER + INCOME + EDUC + AP2, data=apr24.PPP)
apr24.PPP.med.out.y <- mediate(apr24.PPP.med.y, apr24.PPP.out.y, treat = "SORTING", 
                          mediator = "AP2", robustSE = T, sims = 500)

png("figure8.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr24.DEM.med.out.y , main="Democratic Party (April 2024)",cex.main=0.95)
plot(apr24.PPP.med.out.y, main="Conservative Party (April 2024)",cex.main=0.95)
dev.off()

# Appendix: Section A1
# Figure A1
# Originial data downloaded from: https://v-dem.net/data/the-v-dem-dataset/

df <- read.csv("Vdem_appendix.csv")

label_map_01 <- c(
  Equ_Lib = "Individual Liberty",
  Legi_const = "Legislative Constraints\n on the Executive",
  Judi_const = "Judicial Constraints\n on the Executive"
)

df_01 <- df %>%
  dplyr::select(Year, all_of(names(label_map_01))) %>%
  pivot_longer(-Year, names_to = "Variable", values_to = "Value") %>%
  mutate(Variable = label_map_01[Variable])

plot_01 <- ggplot(df_01, aes(x = Year, y = Value, linetype = Variable)) +
  geom_line(size = 1) +
  scale_y_continuous(limits = c(0.5, 1)) +
  scale_x_continuous(limits = c(2000, 2024), breaks = seq(2000, 2024, 4)) +
  labs(
    title = "V-Dem Liberal Democracy Aggregate Indices (0-1)",
    color = NULL
  ) +
  xlab(NULL) + ylab(NULL) +
  theme_minimal(base_family = "Roboto", base_size = 13) +
  theme(
    legend.position = "right",
    legend.text = element_text(size = 9),
    legend.key.width = unit(0.78, "cm"),
    plot.title = element_text(face = "bold", size = 14, hjust = 0.5, margin = margin(b = 10))
  )

print(plot_01)

# Figure A2
label_map_04 <- c(
  academic = "Freedom of Academic\n and Cultural Expression",
  media_censor = "Government Media \nCensorship",
  freedom_exp = "Freedom of Expression\n and Assembly"
)

df_04 <- df %>%
  dplyr::select(Year, academic, media_censor, freedom_exp) %>%
  pivot_longer(-Year, names_to = "Variable", values_to = "Value") %>%
  mutate(Variable = label_map_04[Variable])

plot_04 <- ggplot(df_04, aes(x = Year, y = Value, linetype = Variable)) +
  geom_line(size = 1) +
  scale_y_continuous(limits = c(2, 4)) +
  scale_x_continuous(limits = c(2000, 2024), breaks = seq(2000, 2024, 4)) +
  labs(
    title = "V-Dem Freedom Indices (0-4)",
    x = "",
    y = "",
    color = "") +
  theme_minimal(base_family = "Roboto", base_size = 13) +
  theme(
    legend.position = "right",
    legend.text = element_text(size = 9),
    legend.key.width = unit(0.78, "cm"),
    plot.title = element_text(face = "bold", size = 14, hjust = 0.5, margin = margin(b = 10))
  )

print(plot_04)

# Section A4
# Figure A3
png("figureA3.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1,2))
plot(apr20.DEM.out.y, which = 1, main = "Democratic Party (April 2020) \n Residuals vs. Fitted Plot",
     sub = " ", caption = "")
plot(apr20.MTD.out.y, which = 1, main = "Conservative Party (April 2020) \n Residuals vs. Fitted Plot", 
     sub = " ", caption = "")
dev.off()

# Figure A4
png("figureA4.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1,2))
plot(apr22.DEM.out.y, which = 1, main = "Democratic Party (April 2022) \n Residuals vs. Fitted Plot",
     sub = " ", caption = "")
plot(apr22.PPP.out.y, which = 1, main = "Conservative Party (April 2022) \n Residuals vs. Fitted Plot", 
     sub = " ", caption = "")
dev.off()

# Figure A5
png("figureA5.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1,2))
plot(jun22.DEM.out.y, which = 1, main = "Democratic Party (June 2022) \n Residuals vs. Fitted Plot",
     sub = " ", caption = "")
plot(jun22.PPP.out.y, which = 1, main = "Conservative Party (June 2022) \n Residuals vs. Fitted Plot", 
     sub = " ", caption = "")

# Figure A6
png("figureA6.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1,2))
plot(apr24.DEM.out.y, which = 1, main = "Democratic Party (April 2024) \n Residuals vs. Fitted Plot",
     sub = " ", caption = "")
plot(apr24.PPP.out.y, which = 1, main = "Conservative Party (April 2024) \n Residuals vs. Fitted Plot", 
     sub = " ", caption = "")
dev.off()

# Section A5
# Table A9
summary(apr20.DEM.med.out.y)

# Table A10
summary(apr20.MTD.med.out.y)

# Table A11
summary(apr22.DEM.med.out.y)

# Table A12
summary(apr22.PPP.med.out.y)

# Table A13
summary(jun22.DEM.med.out.y)

# Table A14
summary(jun22.PPP.med.out.y)

# Table A15
summary(apr24.DEM.med.out.y)

# Table A16
summary(apr24.PPP.med.out.y)

# Section A6
# Figure A7
apr20.DEM.sens.out.y.indir <- medsens(apr20.DEM.med.out.y, rho.by = 0.1, effect.type = "indirect", sims = 500)
#summary(apr20.DEM.sens.out.y.indir)
apr20.DEM.sens.out.y.dir <- medsens(apr20.DEM.med.out.y, rho.by = 0.1, effect.type = "direct", sims = 500)
#summary(apr20.DEM.sens.out.y.dir)
par(mfrow = c(1, 2))
plot(apr20.DEM.sens.out.y.indir, sens.par = "rho", ylim = c(-0.15, 0.15))
plot(apr20.DEM.sens.out.y.dir, sens.par = "rho", ylim = c(-0.15, 0.15))

# Figure A8
apr22.DEM.sens.out.y <- medsens(apr22.DEM.med.out.y, rho.by = 0.1, effect.type = "indirect", sims = 500)
# summary(apr22.DEM.sens.out.y)
apr22.DEM.sens.out.y.dir <- medsens(apr22.DEM.med.out.y, rho.by = 0.1, effect.type = "direct", sims = 500)
#summary(apr22.DEM.sens.out.y.dir)
par(mfrow = c(1, 2))
plot(apr22.DEM.sens.out.y, sens.par = "rho", ylim = c(-0.15, 0.15))
plot(apr22.DEM.sens.out.y.dir, sens.par = "rho", ylim = c(-0.15, 0.15))

# Figure A9
jun22.DEM.sens.out.dir <- medsens(jun22.DEM.med.out.y, rho.by = 0.1, effect.type = "direct", sims = 500)
par(mfrow = c(1, 1))
plot(jun22.DEM.sens.out.dir)

# Figure A10
jun22.PPP.sens.out.indir <- medsens(jun22.PPP.med.out.y, rho.by = 0.1, effect.type = "indirect", sims = 500)
par(mfrow = c(1, 1))
plot(jun22.PPP.sens.out.indir, sens.par = "rho")

# Figure A11
apr24.DEM.sens.out.dir <- medsens(apr24.DEM.med.out.y, rho.by = 0.1, effect.type = "direct", sims = 500)
plot(apr24.DEM.sens.out.dir)

# Section A7
# Figure A12
apr24.DEM.ynew <- apr24.DEM %>% filter(!is.na(Ynew))
apr24.DEM.med.ynew <- lm(AP2~SORTING + AGE + GENDER+ INCOME + EDUC, data = apr24.DEM.ynew)
apr24.DEM.out.ynew <- lm(Ynew~SORTING + AGE + GENDER+ INCOME + EDUC + AP2, data = apr24.DEM.ynew)
apr24.DEM.med.out.ynew <- mediate(apr24.DEM.med.ynew , apr24.DEM.out.ynew, 
                             treat = "SORTING", mediator = "AP2", robustSE = F, sims = 500)
#summary(apr24.DEM.med.out.ynew)

apr24.PPPnew <- apr24.PPP %>% filter(!is.na(Ynew))
apr24.PPP.med.ynew <- lm(AP2~SORTING + AGE + GENDER + INCOME + EDUC, data=apr24.PPPnew)
apr24.PPP.out.ynew <- lm(Ynew~SORTING + AGE + GENDER + INCOME + EDUC + AP2, data=apr24.PPPnew)
apr24.PPP.med.out.ynew <- mediate(apr24.PPP.med.ynew, apr24.PPP.out.ynew, treat = "SORTING", 
                             mediator = "AP2", robustSE = TRUE, sims = 1000)
#summary(apr24.PPP.med.out.ynew)

png("figureA12.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr24.DEM.med.out.ynew, main="Democratic Party (April 2024)",cex.main=0.95)
mtext("*“Not sure” responses treated as missing values (NA)", 
      side=1, line=2.5, cex=0.7)
plot(apr24.PPP.med.out.ynew, main="Conservative Party (April 2024)",cex.main=0.95)
dev.off()

# Section A8
# Figure A13
apr20.DEM.med.y10 <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.out.y10 <- lm(Y10 ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.med.out.y10 <- mediate(apr20.DEM.med.y10, apr20.DEM.out.y10, 
                               treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)
apr20.MTD.med.y10 <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.out.y10 <- lm(Y10 ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.med.out.y10 <- mediate(apr20.MTD.med.y10, apr20.MTD.out.y10, 
                               treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

png("figureA13.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr20.DEM.med.out.y10, main = "Democratic Party (April 2020)", cex.main = 0.95)
plot(apr20.MTD.med.out.y10, main = "Conservative Party (April 2020)", cex.main = 0.95)
dev.off()

# Figure A13
apr22.DEM.med.y10 <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.DEM)
apr22.DEM.out.y10 <- lm(Y10 ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.DEM)
apr22.DEM.med.out.y10 <- mediate(apr22.DEM.med.y10, apr22.DEM.out.y10, 
                                 treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)
apr22.PPP.med.y10 <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.PPP)
apr22.PPP.out.y10 <- lm(Y10 ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.PPP)
apr22.PPP.med.out.y10 <- mediate(apr22.PPP.med.y10, apr22.PPP.out.y10, 
                                 treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

png("figureA14.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr22.DEM.med.out.y10, main = "Democratic Party (April 2022)", cex.main = 0.95)
plot(apr22.PPP.med.out.y10, main = "Conservative Party (April 2022)", cex.main = 0.95)
dev.off()

# Figure A15
jun22.DEM.med.y10 <- lm(AP ~ SORTING + AGE + GENDER + INCOME + EDUC, data=jun22.DEM)
jun22.DEM.out.y10 <- lm(Y10 ~ SORTING + AP + AGE + GENDER + INCOME + EDUC, data=jun22.DEM)
jun22.DEM.med.out.y10 <- mediate(jun22.DEM.med.y10, jun22.DEM.out.y10, sims=500, treat = "SORTING", mediator = "AP")
jun22.PPP.med.y10 <- lm(AP ~ SORTING + AGE + GENDER + INCOME + EDUC , data=jun22.PPP)
jun22.PPP.out.y10 <- lm(Y10 ~ SORTING + AP + AGE + GENDER + INCOME + EDUC , data=jun22.PPP)
jun22.PPP.med.out.y10 <- mediate(jun22.PPP.med.y10, jun22.PPP.out.y10, sims=500, treat = "SORTING", mediator = "AP")

png("figureA15.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(jun22.DEM.med.out.y10, main="Democratic Party (June 2022)",cex.main=0.95)
plot(jun22.PPP.med.out.y10, main="Conservative Party (June 2022)",cex.main=0.95)
dev.off()

# Figure A16
apr24.DEM.med.y10 <- lm(AP2~SORTING + AGE + GENDER+ INCOME + EDUC, data=apr24.DEM)
apr24.DEM.out.y10 <- lm(Y10 ~SORTING + AGE + GENDER+ INCOME + EDUC + AP2, data=apr24.DEM)
apr24.DEM.med.out.y10 <- mediate(apr24.DEM.med.y10, apr24.DEM.out.y10, 
                            treat = "SORTING", mediator = "AP2", robustSE = F, sims = 500)
apr24.PPP.med.y10 <- lm(AP2~SORTING + AGE + GENDER + INCOME + EDUC, data=apr24.PPP)
apr24.PPP.out.y10 <- lm(Y10 ~SORTING + AGE + GENDER + INCOME + EDUC + AP2, data=apr24.PPP)
apr24.PPP.med.out.y10 <- mediate(apr24.PPP.med.y10, apr24.PPP.out.y10, treat = "SORTING", 
                            mediator = "AP2", robustSE = T, sims = 500)

png("figureA16.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr24.DEM.med.out.y10, main="Democratic Party (April 2024)",cex.main=0.95)
plot(apr24.PPP.med.out.y10, main="Conservative Party (April 2024)",cex.main=0.95)
dev.off()

# Section A9
# Figure A17
apr20.DEM.med.demodict <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.out.demodict <- lm(DEMODICT ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.med.out.demodict <- mediate(apr20.DEM.med.demodict, apr20.DEM.out.demodict, 
                                 treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

apr20.DEM.med.gtd <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.out.gtd <- lm(GTD ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.med.out.gtd <- mediate(apr20.DEM.med.gtd, apr20.DEM.out.gtd, 
                                      treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

png("figureA17.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr20.DEM.med.out.demodict, main="Democratic Party (April 2020)\nOutcome Variable:\nUndemocratic Attitude (1)", cex.main=0.9)
plot(apr20.DEM.med.out.gtd, main="Democratic Party (April 2020)\nOutcome Variable:\nUndemocratic Attitude (2)", cex.main=0.9)
dev.off()

# Figure A18
apr20.MTD.med.demodict <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.out.demodict <- lm(DEMODICT ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.med.out.demodict <- mediate(apr20.MTD.med.demodict, apr20.MTD.out.demodict, 
                                      treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

apr20.MTD.med.gtd <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.out.gtd <- lm(GTD ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.med.out.gtd <- mediate(apr20.MTD.med.gtd, apr20.MTD.out.gtd, 
                                 treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

png("figureA18.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr20.MTD.med.out.demodict, main="Conservative Party (April 2020)\nOutcome Variable:\nUndemocratic Attitude (1)", cex.main=0.9)
plot(apr20.MTD.med.out.gtd, main="Conservative Party (April 2020)\nOutcome Variable:\nUndemocratic Attitude (2)", cex.main=0.9)
dev.off()

# Figure A19
apr22.DEM.med.demo <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.DEM)
apr22.DEM.out.demo <- lm(DEMODICT ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.DEM)
apr22.DEM.med.out.demo <- mediate(apr22.DEM.med.demo, apr22.DEM.out.demo, 
                                 treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

apr22.DEM.out.gtd <- lm(GTD ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.DEM)
apr22.DEM.med.out.gtd <- mediate(apr22.DEM.med.demo, apr22.DEM.out.gtd, 
                                  treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

apr22.DEM.out.dict <- lm(DICT ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.DEM)
apr22.DEM.med.out.dict <- mediate(apr22.DEM.med.demo, apr22.DEM.out.dict, 
                                 treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

png("figureA19.png",width = 1300, height = 530, res = 210)
par(mfrow = c(1, 3))  
plot(apr22.DEM.med.out.demo, main = "Democratic Party (April 2022)\nOutcome Variable:\nUndemocratic Attitude (1)", cex.main = 0.9)
plot(apr22.DEM.med.out.gtd, main = "Democratic Party (April 2022)\nOutcome Variable:\nUndemocratic Attitude (2)", cex.main = 0.9)
plot(apr22.DEM.med.out.dict, main = "Democratic Party (April 2022)\nOutcome Variable:\nUndemocratic Attitude (3)", cex.main = 0.9)
dev.off()

# Figure A20
apr22.PPP.med.demo <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.PPP)
apr22.PPP.out.demo <- lm(DEMODICT ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.PPP)
apr22.PPP.med.out.demo <- mediate(apr22.PPP.med.demo, apr22.PPP.out.demo, 
                                  treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

apr22.PPP.out.gtd <- lm(GTD ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.PPP)
apr22.PPP.med.out.gtd <- mediate(apr22.PPP.med.demo, apr22.PPP.out.gtd, 
                                 treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

apr22.PPP.out.dict <- lm(DICT ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = apr22.PPP)
apr22.PPP.med.out.dict <- mediate(apr22.PPP.med.demo, apr22.PPP.out.dict, 
                                  treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

png("figureA20.png",width = 1300, height = 530, res = 210)
par(mfrow = c(1, 3))  
plot(apr22.PPP.med.out.demo, main = "Conservative Party (April 2022)\nOutcome Variable:\nUndemocratic Attitude (1)", cex.main = 0.9)
plot(apr22.PPP.med.out.gtd, main = "Conservative Party (April 2022)\nOutcome Variable:\nUndemocratic Attitude (2)", cex.main = 0.9)
plot(apr22.PPP.med.out.dict, main = "Conservative Party (April 2022)\nOutcome Variable:\nUndemocratic Attitude (3)", cex.main = 0.9)
dev.off()

# Figure A21
jun22.DEM.med.demo <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = jun22.DEM)
jun22.DEM.out.demo <- lm(DEMODICT ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = jun22.DEM)
jun22.DEM.med.out.demo <- mediate(jun22.DEM.med.demo, jun22.DEM.out.demo, 
                                  treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

jun22.DEM.out.gtd <- lm(GTD ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = jun22.DEM)
jun22.DEM.med.out.gtd <- mediate(jun22.DEM.med.demo, jun22.DEM.out.gtd, 
                                 treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

jun22.DEM.out.dict <- lm(DICT ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = jun22.DEM)
jun22.DEM.med.out.dict <- mediate(jun22.DEM.med.demo, jun22.DEM.out.dict, 
                                  treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

png("figureA21.png",width = 1300, height = 530, res = 210)
par(mfrow = c(1, 3))  
plot(jun22.DEM.med.out.demo, main = "Democratic Party (June 2022)\nOutcome Variable:\nUndemocratic Attitude (1)", cex.main = 0.9)
plot(jun22.DEM.med.out.gtd, main = "Democratic Party (June 2022)\nOutcome Variable:\nUndemocratic Attitude (2)", cex.main = 0.9)
plot(jun22.DEM.med.out.dict, main = "Democratic Party (June 2022)\nOutcome Variable:\nUndemocratic Attitude (3)", cex.main = 0.9)
dev.off()

# Figure A22
jun22.PPP.med.demo <- lm(AP ~ SORTING + AGE + GENDER + EDUC + INCOME , data = jun22.PPP)
jun22.PPP.out.demo <- lm(DEMODICT ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = jun22.PPP)
jun22.PPP.med.out.demo <- mediate(jun22.PPP.med.demo, jun22.PPP.out.demo, 
                                  treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

jun22.PPP.out.gtd <- lm(GTD ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = jun22.PPP)
jun22.PPP.med.out.gtd <- mediate(jun22.PPP.med.demo, jun22.PPP.out.gtd, 
                                 treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

jun22.PPP.out.dict <- lm(DICT ~ AP + SORTING + AGE + GENDER + EDUC + INCOME , data = jun22.PPP)
jun22.PPP.med.out.dict <- mediate(jun22.PPP.med.demo, jun22.PPP.out.dict, 
                                  treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

png("figureA22.png",width = 1300, height = 530, res = 210)
par(mfrow = c(1, 3))  
plot(jun22.PPP.med.out.demo, main = "Conservative Party (June 2022)\nOutcome Variable:\nUndemocratic Attitude (1)", cex.main = 0.9)
plot(jun22.PPP.med.out.gtd, main = "Conservative Party (June 2022)\nOutcome Variable:\nUndemocratic Attitude (2)", cex.main = 0.9)
plot(jun22.PPP.med.out.dict, main = "Conservative Party (June 2022)\nOutcome Variable:\nUndemocratic Attitude (3)", cex.main = 0.9)
dev.off()

# Figure A23
apr24.DEM.med.demo <- lm(AP2 ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr24.DEM)
apr24.DEM.out.demo <- lm(DEMODICT ~ AP2 + SORTING + AGE + GENDER + EDUC + INCOME , data = apr24.DEM)
apr24.DEM.med.out.demo <- mediate(apr24.DEM.med.demo, apr24.DEM.out.demo, 
                                  treat = "SORTING", mediator = "AP2", robustSE = T, sims = 500)

apr24.DEM.out.gtd <- lm(GTD ~ AP2 + SORTING + AGE + GENDER + EDUC + INCOME , data = apr24.DEM)
apr24.DEM.med.out.gtd <- mediate(apr24.DEM.med.demo, apr24.DEM.out.gtd, 
                                 treat = "SORTING", mediator = "AP2", robustSE = T, sims = 500)

apr24.DEM.out.dict <- lm(DICT ~ AP2 + SORTING + AGE + GENDER + EDUC + INCOME , data = apr24.DEM)
apr24.DEM.med.out.dict <- mediate(apr24.DEM.med.demo, apr24.DEM.out.dict, 
                                  treat = "SORTING", mediator = "AP2", robustSE = T, sims = 500)

png("figureA23.png",width = 1300, height = 530, res = 210)
par(mfrow = c(1, 3))  
plot(apr24.DEM.med.out.demo, main = "Democratic Party (April 2024)\nOutcome Variable:\nUndemocratic Attitude (1)", cex.main = 0.9)
plot(apr24.DEM.med.out.gtd, main = "Democratic Party (April 2024)\nOutcome Variable:\nUndemocratic Attitude (2)", cex.main = 0.9)
plot(apr24.DEM.med.out.dict, main = "Democratic Party (April 2024)\nOutcome Variable:\nUndemocratic Attitude (3)", cex.main = 0.9)
dev.off()

# Figure A24
apr24.PPP.med.demo <- lm(AP2 ~ SORTING + AGE + GENDER + EDUC + INCOME , data = apr24.PPP)
apr24.PPP.out.demo <- lm(DEMODICT ~ AP2 + SORTING + AGE + GENDER + EDUC + INCOME , data = apr24.PPP)
apr24.PPP.med.out.demo <- mediate(apr24.PPP.med.demo, apr24.PPP.out.demo, 
                                  treat = "SORTING", mediator = "AP2", robustSE = T, sims = 500)

apr24.PPP.out.gtd <- lm(GTD ~ AP2 + SORTING + AGE + GENDER + EDUC + INCOME , data = apr24.PPP)
apr24.PPP.med.out.gtd <- mediate(apr24.PPP.med.demo, apr24.PPP.out.gtd, 
                                 treat = "SORTING", mediator = "AP2", robustSE = T, sims = 500)

apr24.PPP.out.dict <- lm(DICT ~ AP2 + SORTING + AGE + GENDER + EDUC + INCOME , data = apr24.PPP)
apr24.PPP.med.out.dict <- mediate(apr24.PPP.med.demo, apr24.PPP.out.dict, 
                                  treat = "SORTING", mediator = "AP2", robustSE = T, sims = 500)

png("figureA24.png",width = 1300, height = 530, res = 210)
par(mfrow = c(1, 3))  
plot(apr24.PPP.med.out.demo, main = "Conservative Party (April 2024)\nOutcome Variable:\nUndemocratic Attitude (1)", cex.main = 0.9)
plot(apr24.PPP.med.out.gtd, main = "Conservative Party (April 2024)\nOutcome Variable:\nUndemocratic Attitude (2)", cex.main = 0.9)
plot(apr24.PPP.med.out.dict, main = "Conservative Party (April 2024)\nOutcome Variable:\nUndemocratic Attitude (3)", cex.main = 0.9)
dev.off()

# Section A10
# Figure A25
apr20.DEM.med.s2 <- lm(AP ~ SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.out.s2 <- lm(Y ~ AP + SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.med.out.s2 <- mediate(apr20.DEM.med.s2, apr20.DEM.out.s2, 
                              treat = "SORTING2", mediator = "AP", robustSE = TRUE, sims = 500)
apr20.MTD.med.s2 <- lm(AP ~ SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.out.s2 <- lm(Y ~ AP + SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.med.out.s2 <- mediate(apr20.MTD.med.s2, apr20.MTD.out.s2, 
                              treat = "SORTING2", mediator = "AP", robustSE = T, sims = 500)

png("figureA25.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr20.DEM.med.out.s2, main="Democratic Party (April 2020)",cex.main=0.95)
mtext("*Partisan Sorting with different Political Ideology classification", 
      side=1, line=2.5, cex=0.7)
mtext("(Liberal: 0~3; Independent: 4~6; Conservative: 7~10)", 
      side=1, line=3, cex=0.7)
plot(apr20.MTD.med.out.s2, main="Conservative Party (April 2020)", cex.main=0.95)
dev.off()

# Figure A26
apr22.DEM.med.s2 <- lm(AP ~ SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr22.DEM)
apr22.DEM.out.s2 <- lm(Y ~ AP + SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr22.DEM)
apr22.DEM.med.out.s2 <- mediate(apr22.DEM.med.s2, apr22.DEM.out.s2, 
                                treat = "SORTING2", mediator = "AP", robustSE = TRUE, sims = 500)
apr22.PPP.med.s2 <- lm(AP ~ SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr22.PPP)
apr22.PPP.out.s2 <- lm(Y ~ AP + SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr22.PPP)
apr22.PPP.med.out.s2 <- mediate(apr22.PPP.med.s2, apr22.PPP.out.s2, 
                               treat = "SORTING2", mediator = "AP", robustSE = T, sims = 500)

png("figureA26.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr22.DEM.med.out.s2, main="Democratic Party (April 2022)",cex.main=0.95)
mtext("*Partisan Sorting with different Political Ideology classification", 
      side=1, line=2.5, cex=0.7)
mtext("(Liberal: 0~3; Independent: 4~6; Conservative: 7~10)", 
      side=1, line=3, cex=0.7)
plot(apr22.PPP.med.out.s2, main="Conservative Party (April 2022)", cex.main=0.95)
dev.off()

# Figure A27
jun22.DEM.med.s2 <- lm(AP ~ SORTING2 + AGE + GENDER + EDUC + INCOME , data = jun22.DEM)
jun22.DEM.out.s2 <- lm(Y ~ AP + SORTING2 + AGE + GENDER + EDUC + INCOME , data = jun22.DEM)
jun22.DEM.med.out.s2 <- mediate(jun22.DEM.med.s2, jun22.DEM.out.s2, 
                                treat = "SORTING2", mediator = "AP", robustSE = TRUE, sims = 500)
jun22.PPP.med.s2 <- lm(AP ~ SORTING2 + AGE + GENDER + EDUC + INCOME , data = jun22.PPP)
jun22.PPP.out.s2 <- lm(Y ~ AP + SORTING2 + AGE + GENDER + EDUC + INCOME , data = jun22.PPP)
jun22.PPP.med.out.s2 <- mediate(jun22.PPP.med.s2, jun22.PPP.out.s2, 
                                treat = "SORTING2", mediator = "AP", robustSE = T, sims = 500)

png("figureA27.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(jun22.DEM.med.out.s2, main="Democratic Party (June 2022)",cex.main=0.95)
mtext("*Partisan Sorting with different Political Ideology classification", 
      side=1, line=2.5, cex=0.7)
mtext("(Liberal: 0~3; Independent: 4~6; Conservative: 7~10)", 
      side=1, line=3, cex=0.7)
plot(jun22.PPP.med.out.s2, main="Conservative Party (June 2022)", cex.main=0.95)
dev.off()

# Figure A28
apr24.DEM.med.s2 <- lm(AP2 ~ SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr24.DEM)
apr24.DEM.out.s2 <- lm(Y ~ AP2 + SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr24.DEM)
apr24.DEM.med.out.s2 <- mediate(apr24.DEM.med.s2, apr24.DEM.out.s2, 
                                treat = "SORTING2", mediator = "AP2", robustSE = TRUE, sims = 500)
apr24.PPP.med.s2 <- lm(AP2 ~ SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr24.PPP)
apr24.PPP.out.s2 <- lm(Y ~ AP2 + SORTING2 + AGE + GENDER + EDUC + INCOME , data = apr24.PPP)
apr24.PPP.med.out.s2 <- mediate(apr24.PPP.med.s2, apr24.PPP.out.s2, 
                                treat = "SORTING2", mediator = "AP2", robustSE = T, sims = 500)

png("figureA28.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr24.DEM.med.out.s2, main="Democratic Party (April 2024)",cex.main=0.95)
mtext("*Partisan Sorting with different Political Ideology classification", 
      side=1, line=2.5, cex=0.7)
mtext("(Liberal: 0~3; Independent: 4~6; Conservative: 7~10)", 
      side=1, line=3, cex=0.7)
plot(apr24.PPP.med.out.s2, main="Conservative Party (April 2024)", cex.main=0.95)
dev.off()

# Section A11
# Figure A29

apr20.DEM.med.long <- lm(AP ~ SORTINGlong + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.out.long <- lm(Y ~ AP + SORTINGlong + AGE + GENDER + EDUC + INCOME , data = apr20.DEM)
apr20.DEM.med.out.long <- mediate(apr20.DEM.med.long, apr20.DEM.out.long, 
                                  treat = "SORTINGlong", mediator = "AP", robustSE = TRUE, sims = 500)
apr20.MTD.med.long <- lm(AP ~ SORTINGlong + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.out.long <- lm(Y ~ AP + SORTINGlong + AGE + GENDER + EDUC + INCOME , data = apr20.MTD)
apr20.MTD.med.out.long <- mediate(apr20.MTD.med.long, apr20.MTD.out.long, 
                                  treat = "SORTINGlong", mediator = "AP", robustSE = T, sims = 500)

png("figureA29.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr20.DEM.med.out.long, main="Democratic Party (April 2020)",cex.main=0.95)
mtext("*Including significant third apr20ties (Justice Party and People Party)", 
      side=1, line=2.5, cex=0.7)
plot(apr20.MTD.med.out.long, main="Conservative Party (April 2020)", cex.main=0.95)
dev.off()

# Section A12
# Figure A30
apr24.DEM.med.ap1 <- lm(AP~SORTING + AGE + GENDER+ INCOME + EDUC, data=apr24.DEM)
apr24.DEM.out.ap1 <- lm(Y~SORTING + AGE + GENDER+ INCOME + EDUC + AP, data=apr24.DEM)
apr24.DEM.med.out.ap1 <- mediate(apr24.DEM.med.ap1, apr24.DEM.out.ap1, 
                          treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

apr24.PPP.med.ap1 <- lm(AP~SORTING + AGE + GENDER+ INCOME + EDUC, data=apr24.PPP)
apr24.PPP.out.ap1 <- lm(Y~SORTING + AGE + GENDER+ INCOME + EDUC + AP, data=apr24.PPP)
apr24.PPP.med.out.ap1 <- mediate(apr24.PPP.med.ap1, apr24.PPP.out.ap1, 
                                 treat = "SORTING", mediator = "AP", robustSE = T, sims = 500)

png("figureA30.png", width = 1050, height = 600, res = 150)
par(mfrow = c(1, 2))
plot(apr24.DEM.med.out.ap1, main="Democratic Party (April 2024)",cex.main=0.95)
mtext("*Affective polarization as the difference in average feeling thermometer", 
      side=1, line=2.5, cex=0.7)
plot(apr24.PPP.med.out.ap1, main="Conservative Party (April 2024)",cex.main=0.95)
mtext("ratings of major Democratic/Conservative Party figures", 
      side=1, line=2.5, cex=0.7)
dev.off()

# Section A13
# Figure A31
jun22.PPP.sort <- jun22.PPP %>% filter(SORTING==1)
apr24.PPP.sort <- apr24.PPP %>% filter(SORTING==1)

mean1 <- mean(jun22.PPP.sort$Y)
mean2 <- mean(apr24.PPP.sort$Y)
se1 <- sd(jun22.PPP.sort$Y) / sqrt(length(jun22.PPP.sort$Y))
se2 <- sd(apr24.PPP.sort$Y) / sqrt(length(apr24.PPP.sort$Y))

df_means <- data.frame(
  group = c("June 2022", "April 2024"),
  mean = c(mean1, mean2),
  se = c(se1, se2)
)
df_means$group <- factor(df_means$group, levels = c("June 2022", "April 2024"))

ggplot(df_means, aes(x = group, y = mean, fill = group)) +
  geom_bar(stat = "identity", width = 0.5, alpha = 0.7) +
  geom_errorbar(aes(ymin = mean - 1.96 * se, ymax = mean + 1.96 * se),
                width = 0.2) +
  theme_minimal() +
  labs(title = "Comparison of Average Undemocratic Attitudes\namong Sorted Conservative Party Supporters",
       x = "", y = "Mean of Average Undemocratic Attitudes") +
  theme(
    plot.title = element_text(face = "bold",hjust = 0.5),
    legend.position = "none"
  )

ggsave("figureA31.jpeg",
       width = 6, height = 5, dpi = 300)

# Figure A32
apr20.MTD.nots <- apr20.MTD %>% filter(SORTING==0)
apr22.PPP.nots <- apr22.PPP %>% filter(SORTING==0)
jun22.PPP.nots <- jun22.PPP %>% filter(SORTING==0)
apr24.PPP.nots <- apr24.PPP %>% filter(SORTING==0)

y1 <- apr20.MTD.nots$Y
y2 <- apr22.PPP.nots$Y
y3 <- jun22.PPP.nots$Y
y4new <- apr24.PPP.nots$Ynew

g1 <- rep("April_2020", length(y1))
g2 <- rep("April_2022", length(y2))
g3 <- rep("June_2022", length(y3))
g4new <- rep("April_2024", length(y4new))

df_anova_new <- data.frame(
  Y = c(y1, y2, y3, y4new),
  group = factor(c(g1, g2, g3, g4new),
                 levels = c("April_2020", "April_2022", "June_2022", "April_2024"))
)

df_summary_new <- df_anova_new %>%
  group_by(group) %>%
  summarise(
    mean_Y = mean(Y, na.rm=T),
    se_Y = sd(Y, na.rm=T) / sqrt(n())
  )

summary(aov(Y ~ group, data = df_anova_new))

ggplot(df_summary_new, aes(x = group, y = mean_Y, fill = group)) +
  geom_bar(stat = "identity", width = 0.6, alpha = 0.8) +
  geom_errorbar(aes(ymin = mean_Y - 1.96 * se_Y,
                    ymax = mean_Y + 1.96 * se_Y),
                width = 0.2) +
  theme_minimal() +
  labs(
    title = "Comparison of Average Undemocratic Attitudes\n among Unsorted Conservative Party Supporters",
    x = "",
    y = "Mean of Average Undemocratic Attitudes"
  ) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    legend.position = "none"
  )

ggsave("figureA32.jpeg",
       width = 5.5, height = 4.8, dpi = 300)
